#import libraries
library(irr)
library(psych)
library(dplyr)
library(car)
library(effectsize)
library(emmeans)

#load data
data <- read.csv('Study 4a.csv')
View(data)

#recode categorical vars with values
data$gender_coded <- recode_factor(data$gender, 
                                   '1' = 'male', 
                                   '2' = 'female',
                                   '4' = 'other',
                                   .default = 'missing')


data$exp_specific_coded <- recode_factor(data$exp_specific,
                                         '1' = 'experienced',
                                         '2' = 'inexperienced')

data$exp_general_coded <- recode_factor(data$exp_general,
                                        '1' = 'experienced',
                                        '2' = 'inexperienced')


#get the sample information
##condition
condition_count<-table(data$condition) #low-inequality|low-income = 118 vs. low-inequality|middle-income = 122 vs. high-inequality = 119
condition_count

##gender
gender_count<-table(data$gender_coded) #male = 182 vs. female = 173, 1 prefer not to say, 3 missing
gender_count

##age
summary(data$age) #36.63
sd(data$age) #12.39

##exp
table(data$exp_specific_coded) #w/exp = 21 vs. w/o exp = 338
table(data$exp_general_coded) #w/exp = 58 vs. w/o exp = 299


#manipulation check
##calculate the manipulation check of inequality 
data <- data %>%
  mutate(manipulation_check = (check_1 + check_2)/2)

##reliability 
inequality_alpha <- data[, c('check_1', 'check_2')]
cor(inequality_alpha) # correlation = .69

##check
aggregate(data$manipulation_check ~ data$condition, FUN = mean) 
aggregate(data$manipulation_check ~ data$condition, FUN = sd)
##low-inequality|low-income = 3.32(1.16) vs. low-inequality|middle-income = 3.05(1.23) vs. high-inequality = 5.21(1.22)


aov_check <- aov(manipulation_check ~ condition, data = data)
summary(aov_check)  #F(2,356) = 114.2, p < .001


#reliability for coding
##coding for SES inference
kappa2(data[,c('SES_inference_coder1','SES_inference_coder2')]) #Cohen's Kappa = .66
##addtional test
data_test <- data[!is.na(data$SES_inference) & data$SES_inference != 0, ]
cor_test <- cor.test(data_test$SES, data_test$SES_inference,use = "complete.obs") #p = .75, p < .001
cor_test

##conding for trustworthiness
icc(data[,c('trustworthiness_coder1','trustworthiness_coder2')], 
    model = "twoway", type = "consistency", unit = "average") #ICC-consistency = .766
icc(data[,c('trustworthiness_coder1','trustworthiness_coder2')],
    model = "twoway", type = "agreement", unit = "average") #ICC-agreement = .735


#calculate the mean of the final scores for coding
##final score for SES inference
data <- data %>%
  mutate(SES_inference = (SES_inference_coder1 + SES_inference_coder2)/2)

##final score for trustworthiness
data <- data %>%
  mutate(trustworthiness = (trustworthiness_coder1 + trustworthiness_coder2)/2)


#main effect of inequality on willingness
##descriptive
aggregate(data$willingness ~ data$condition, FUN = mean)
aggregate(data$willingness ~ data$condition, FUN = sd)
##low-inequality|low-income = 4.38(1.74) vs. low-inequality|middle-income = 5.39(1.18) vs. high-inequality = 5.30(1.56)

#ANOVA
aov_willingness <- aov(willingness ~ condition, data)
summary(aov_willingness) #F(2, 356) = 32.16, p < .001
effectsize(aov_willingness, partial = TRUE) #eta-sq = .15
emm_willingness <- (emmeans(aov_willingness, ~ condition))
pairs(emm_willingness, adjust = 'None')
eff_size(emm_willingness, sigma = sigma(aov_willingness), edf = df.residual(aov_willingness))

#main effect of inequality on perceived trustworthiness (5-point scale)
##descriptive
aggregate(data$trustworthiness ~ data$condition, FUN = mean)
aggregate(data$trustworthiness ~ data$condition, FUN = sd)
##low-inequality|low-income = 2.75(.94) vs. low-inequality|middle-income = 4.07(1.03) vs. high-inequality = 3.53(.95)

##ANOVA
aov_trust <- aov(trustworthiness ~ condition, data)
summary(aov_trust) #F(2, 356) = 54.73, p < .001
effectsize(aov_trust, partial = TRUE) #eta-sq = .24
emm_trust <- (emmeans(aov_trust, ~ condition))
pairs(emm_trust, adjust = 'None')
eff_size(emm, sigma = sigma(aov_trust), edf = df.residual(aov_trust))


#mediation analysis
process(data = data, 
        y = "willingness", 
        x = "condition_code", 
        m = c("trustworthiness"),
        model = 4,
        total = 1,
        boot = 5000,
        modelbt = 1,
        seed = 123456)

##main effect of inequality on SES inference
##descriptive
aggregate(data$SES_inference ~ data$condition, FUN = mean) #
aggregate(data$SES_inference ~ data$condition, FUN = sd)
##low-inequality|low-income = -.25(.45) vs. low-inequality|middle-income = .46(.46) vs. high-inequality = .0042(.45)

##ANOVA
aov_SES <- aov(SES_inference ~ condition, data)
summary(aov_SES) #F(2, 356) = 78.41, p < .001
effectsize(aov_SES, partial = TRUE) #eta-sq = .31
emm_SES <- (emmeans(aov_SES, ~ condition))
pairs(emm_SES, adjust = 'None')
eff_size(emm_SES, sigma = sigma(aov_SES), edf = df.residual(aov_SES))

##serial mediation
process(data = data, 
        y = "willingness", 
        x = "condition_code", 
        m = c("SES_inference","trustworthiness"),
        model = 6,
        total = 1,
        boot = 5000,
        modelbt = 1,
        seed = 123456)


##simple mediation
process(data = data, 
        y = "willingness", 
        x = "condition_code", 
        m = c("trustworthiness"),
        model = 4,
        total = 1,
        boot = 5000,
        modelbt = 1,
        seed = 123456)